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ABSTRACT 

Some new iterative methods for numerical solution of mixed finite element 
approximation of Stokes problem are presented. The idea is the use of proper 
preconditioning for the conjugate gradient algorithm. A particular case gives 
a variant of the Arrow -Hurwicz method. 
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1. STATEMENT OF THE PROBLEM 

Let us consider a polygonal domain Q. aR n (n=2 or 3) of regular boundary dQ. — T. 
Let us denote V = [v E (Hl(0)) n }, Q = L 2 (H), e{v) = (e y (v)) and 

The Stokes problem for fluid flow is 
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j - v ZU£-j e vM + (Vp)t = f u 1 < i < n in ft, 

I u E V, p e Q, 

Where / = (f lt ...,/ n ) is a given function inQ n and v the viscosity of the fluid. We 
look for a velocity vector u and a static pressure p. The problem (1.1) - (1.2) can be 
formulated as the saddle point problem. 

Inf VEV Sup qEQ L{v,q) = Inf VEV Sup qEQ (J(v) - (q.V.v)) (1.3) 

Wherey(i/) = \v J n e(v): e(v)dx - f Q f.vdx = \a{v,v)-(f,v) 

With the condition V. v = 0 we obtain 

cl(v, v) = v j a e(v): e(v)dx — v j a Vv.Vv dx (1.4) 

Taking the equilibrium condition for (1.3) gives: 

{v f a Vu: Vv dx - f a f. v dx - / fl pV. v dx = 0 V v G V, 
f a V.uq dx = 0 Vq G Q. 

This is clearly a variational formulation of (1.1) - (1.2). Let V h czV and Qh^Q be 
two families of finite element spaces approximations V and Q respectively, A and B 
the discrete operators approaching (— vA) and V respectively in these approximation. 
The discretized Stokes problem is to find (u, p) G V h x Q h , solution of the following 
symmetric and indefinite system 



(Au + B t p=fina (l6) 
{ Bu = 0 in ft 



VB 0-1' Lpl LOT 
The problem (1.6) takes the following form 



(1.7) 



A is indefinite symmetric matrix with a particular structure. The discrete primal 
problem obtained by restricting to the divergence-free subspace can be written as 

Inf veVoh \(A 0 v,v) - (J 0 ,v) (1.8) 



V oh = {veV h \Bv = 0} (1.9) 
Where A 0 and f 0 are respectively the restriction of A and f to V oh , (1.8) is equivalent 
to the minimization problem. 
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The discrete dual problem is 

Inf qeQh \CDq,q)-(r,q) (1.10) 

Where D = BA~ x B l and /* = BA~ X f, this problem is also a standard quadratic 
problem although D = BA~ 1 B t cannot be computed explicitly as A' 1 is a full matrix. 
We shall develop algorithms taking into account the special structure of A to obtain a 
general family of iterative methods some of which will be explicited and analysed. 
Before doing so we recall the application of classical optimization techniques to the 
primal and dual problems (1.8) and (1.10). 

2. CLASSICAL SOLUTION METHODS 

2. 1. Primal Problem 

The Navier-Stokes and Stokes problem require a preconditioner yielding a 
divergence-free solution and a divergence-free descent direction at each iteration [2, 
7]. An efficient example of this preconditioner can be built throught the block 
relaxation method, the advantage of which is to solve a discretized problem in a small 
subregion. We have described in previous paper (cf [1], [4], [5]) such method. Let us 
denote S this preconditioning operator, the P.C.G. algorithm becomes. 

2.1.1. Alg. 2.1: 

Step 1: Select an initial divergence-free solution u° 
Step 2: Compute the divergence-free descent direction 

z n = 5" V (2.2) 
Where g n = Au n - f, and compute 

<D" = z n + Bn® 71 ' 1 so that A® 71 !® 11 ' 1 (2.3) 

Step 3: Compute a n ,u n+1 and g n+1 by 

u n+1 = u"-a n <t>" 
g n+1 = Au n+1 - f 
a n defined by the condition g n+1 ±<$> n . 

This is in fact the usual P.C.G method on the divergence-free subspace. 

2.2. Dual Problem 

We remember in the following the C.G Uzawa algorithm for resolution of the dual 
problem. 

2.2.1. Alg. 2.2: 

Step 1: Let p° G Q h ,u° = A' 1 ^ - 5 f p 0 ), suppose u n known. 
Step 2: Compute the descent direction 

z n = [z^.Zp) = (-A~ 1 B t Bu n ,Bu n ) and compute 
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Step 3: Compute a n ,u n+1 and p n+1 by 



(2.5) 



(2.6) 



(u n+1 = u n - a n <$>1 

[ p n+l = p n _ an0 n (2-7) 

This algorithm requires to solve exactly the linear system An = f — B t p at each 
iteration. If we solve approximatively this system, we obtain Arrow-Hurwicz's 
algorithm to find the saddle point. The steepest descent method is obtained for B n = 
0; for a = 0 constant we obtain the Uzawa algorithm. 

3. GENERAL FORMULATION 

The principal idea of this method is to combine the two algorithms Alg. 2. 1 and alg. 
2.2 for solving mutually the primal and the dual problems. The system (1.7) is 
indefinite, the standard C.G method yields a divergent iterative method. However 
with a good preconditioner and proper descent direction we can obtain a convergent 
iterative method, which coincides with a variant of Arrow-Hurwicz algorithm (cf. 
[3]). The next algorithm is interesting when the projection on a divergence-free 
subspace is difficult or very expensive with the preconditioner used. Let S be a 
preconditioning operator of A (cf. [1]), R n and R n are the residuals respectively 
defined by 



where r u n = Au n + 5 f p n - / (3.1) 
and 

R n = [ r n] where r p n = Bu n (3.2) 
The step descent directions are defined as solutions of the following problems 

Sz n = R n (3.3) 

Sz n = R n (3.4) 
Those directions are defined to minimize respectively the resuduals of the primal 
and the dual problems. 

3.1. G. "Primal-Dual" Algorithm 

Step 1: Select an initial solution(u°,p°). 

http://www.iaeme.com/IJARET/index.asp (^^^) editor@iaeme.com 



Ndzana Benoit, Biya Motto And Lekini Nkodo Claude Bernard 
Step 2: Solve Sz 11 = R n 

a n is defined by a condition in order to minimize the residual R n , 



r u 2 =u n -a n z£ (3.5) 
Compute u~ and p~ 



u 2 = u' — a n z£ 

p — = p n_ anZ n 



r u 2 = Au 2 + B t p z - / 
r p ~ = Bu^r 
Step 3: Solve Sz n = Rn+i 

a n defined by a condition in order to minimize the residualftn+i, i.e. 

2 

Compute u n+1 and p n+1 



= u 2 - a n Zu 
P n+1 = P~ ~ a n Zp 



r r n+l _ ^n+l + 5 t p n+l . 

( r p n+1 = Bu n+1 



(3.6) 



(3.7) 



(3.8) 



(3.9) 



(3.10) 



3.2. Primal-Dual Algorithm 

Since every one of the steps 2 and 3 correspond to the gradient algorithm for 
minimization of the residual for a quadratic problem, we can introduce the 
orthogonality relations associated to the primal problem matrix A, and the 
approximate dual problem matrix D = BS~ 1 B t , we obtain the following algorithm. 

3.2.1. Alg. 3.1. 

Step 1: Select an initial solution (u°,p°) 
Step 2: Compute B n and O n by a condition 
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r u 2 ±02 -1 

0 n = z n + p n $>n-l 



Compute a n , U n+i and pn+i ; a n defined by a condition 



= p n - a n z£ 



--Am + B*p : 
r n 2 = Bu^~ 



(3-11) 



(3.12) 
(3.13) 



(3.14) 



Step 3: Solve Sz n = Rn+i 
Compute B n and O n by a condition 

r n + l ±0 n-l (315) 

<D n = z n + ^ n <!)„_! 

Compute a n , u n+1 , p n+1 . a n defined by a condition 

r" +1 ±0£ (3.16) 



r u n+1 = Au n+1 + B l p n+1 - f 
r p n+1 = Bu n+1 

Doing some iterations of Step 2 before moving to Step 3, we obtain a convergence 
of the primal variable corresponding to a minimum in v of the lagrangien L(y, q) and 
we pull back on Uzawa's method (Alg. 2.2). 

On the contrary if we obtain convergence of the approximate dual variable and a 
divergence-free solution after some iteration of Step 3, the whole process can be 
reduced to the primal algorithm (Alg. 2.1). With a good preconditioner we can obtain 
easily the divergence-free condition and we find once again the alg. 2.1. If the 
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preconditioner yields rapidly the convergence of the primal problem we find once 
again the Uzawa algorithm. The study of the convergence depends of the 
preconditioner (cf. R. Aboulaich [1]). The convergence is illustrated in Figure 1. In 
the following we present a particular example of preconditioning operator and we find 
a variant of Arrow-Hurwicz algorithm. 

3.3. A Particular Case of Preconditioning 



Let us denote S = [jj ^] , R n = [ r g] and R n = [ r ° n ] 



The G. "P-D" algorithm becomes: 
Alg. 3.2 

Step 1: The same 
Step 2: Solve 




(3.19) 



And we compute 




(3.20) 




(3.21) 



Step 3: Solve 




(3.22) 



And we compute 



(3.23) 



|u n+1 = u 2 - a n z" 



(3.24) 



[p n+1 = V 2 ~ a nZp 



We obtain 



u 2 = u n - a n S 1 r£ 



(3.25) 



u 1 



(3.26) 
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P n+1 = P n ~ a n r£ (3.27) 



We find a variant of Arrow-Hurwicz algorithm. A must practical variante of the 
previous algorithm is obtained to compute parallel the two descente directions z and z, 
the algorithm becomes: 

u n+l =u n_ an 5-V M n - OnS-WBll 11 , 

p n + l =p n + anBy n C 3 - 2 *) 

Also we can use the direction g n = z n + yz n where 

L 0 J L -Bu n J 

The iteration (4.12) becomes 



(3.29) 



(pn+l _ pti _ CCn B u n i 

cc n defined to minimize the residual R n = ^J, so that R n+1 lg n , and we have the 

possibility to use the P.C.G algorithm for the saddle point problem, with the recherche 
direction ifj n = g n + ^ n ^ n_1 such that (Aifj n+1 ,ifj n ) = 0. The problem is to define 
properly the parameter y, when y is constant we obtain convergence of the algorithm 
but is difficult to compute the optimal parameter. 
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